SUBROUTINE DG_CONT
!-----------------------------------------------------------------------
!
! Definition von DG:  Kontinuitaetsgleichung
!

      use physco,   only : z1, z2, pi
      use primvar,  only : DG, X, XZ, XA, MD, MR, MUr, &
                           DG_i, DG_im1, DG_ip1, DG_im2, DG_ip2
      use config,   only : np, rkind
      use global,   only : tst, zz, relax
      use advecvar, only : rho_adv, drho_adv_dim2, drho_adv_dim1, drho_adv_di, drho_adv_dip1
      use geomvar,  only : S_flux, S_vol, S_volA


      implicit none

      integer :: i


!-----------------------------------------------------------------------
!    Randbedingungen
!-----------------------------------------------------------------------

      ! innere Pseudozellen: i=1 & i=2
!~       if ( X(MUr,3) * tst / X(MR,3) > +1.e-30_rkind ) then  ! Einstroemung ueber Innenrand

         DG(MD,MD,DG_i,1)      =  z1
         DG(MD,MD,DG_i,2)      =  z1

!~       else
!~
!~          DG(MD,MD,DG_i,1)      =  z1
!~          DG(MD,MD,DG_ip1,1)    = -z1
!~          DG(MD,MD,DG_i,2)      =  z1
!~          DG(MD,MD,DG_ip1,2)    = -z1
!~
!~       end if


      ! aeussere Pseudozellen + ueberzaehlige skalare Zelle: i=np, i=np-1, i=np-2
!~       if ( X(MUr,np-2) * tst / X(MR,np-2) < -1.e-30_rkind ) then  ! Einstroemung ueber Aussenrand

         DG(MD,MD,DG_i,np)     =  z1
         DG(MD,MD,DG_i,np-1)   =  z1
         DG(MD,MD,DG_i,np-2)   =  z1

!~       else
!~
!~          DG(MD,MD,DG_i,np)     =  z1
!~          DG(MD,MD,DG_im1,np)   = -z1
!~          DG(MD,MD,DG_i,np-1)   =  z1
!~          DG(MD,MD,DG_im1,np-1) = -z1
!~          DG(MD,MD,DG_i,np-2)   =  z1
!~          DG(MD,MD,DG_im1,np-2) = -z1
!~
!~       end if


!-----------------------------------------------------------------------
!    Restlicher Bereich
!-----------------------------------------------------------------------

      do i=3,np-3

DG(MD,MR,DG_i,i)=-(pi*X(MD,i)*(X(MR,i) + X(MR,i+1))) + pi*X(MD,i)*(X(MR,i+1) - X(MR,i)) - (-(pi*(XA(MR,i) + &
X(MR,i))) - pi*(X(MR,i) - XA(MR,i)) + pi*tst*XZ(MUr,i)*z2 *zz )*rho_adv(i)
DG(MD,MR,DG_ip1,i)=pi*X(MD,i)*(X(MR,i) + X(MR,i+1)) + pi*X(MD,i)*(X(MR,i+1) - X(MR,i)) + (-(pi*(XA(MR,i+1) + &
X(MR,i+1))) - pi*(X(MR,i+1) - XA(MR,i+1)) + pi*tst*XZ(MUr,i+1)*z2 *zz )*rho_adv(i+1)
DG(MD,MD,DG_im2,i)=-(S_flux(i)*zz *drho_adv_dim2(i))
DG(MD,MD,DG_im1,i)=-(S_flux(i)*zz *drho_adv_dim1(i)) + S_flux(i+1)*zz *drho_adv_dim2(i+1)
DG(MD,MD,DG_i,i)=S_vol(i) - S_flux(i)*zz *drho_adv_di(i) + S_flux(i+1)*zz *drho_adv_dim1(i+1)
DG(MD,MD,DG_ip1,i)=-(S_flux(i)*zz *drho_adv_dip1(i)) + S_flux(i+1)*zz *drho_adv_di(i+1)
DG(MD,MD,DG_ip2,i)=S_flux(i+1)*zz *drho_adv_dip1(i+1)
DG(MD,MUr,DG_i,i)=-(pi*XZ(MR,i)*tst*z2 *zz *rho_adv(i))
DG(MD,MUr,DG_ip1,i)=pi*XZ(MR,i+1)*tst*z2 *zz *rho_adv(i+1)

      end do


END SUBROUTINE DG_CONT
